Example Usage

%matplotlib inline
%config InlineBackend.figure_formats = ['svg']
import chalc as ch
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams["animation.html"] = "jshtml"
import h5py, warnings

For our data we sample 100 points on a circle with some noise and 100 points from inside the unit disk.

rng = np.random.default_rng(40)
num_points_circle = 100
num_points_disk = 100
mean = [0, 0]
cov = np.eye(2)*0.01
circle = np.array([[np.sin(2*np.pi*t), np.cos(2*np.pi*t)] for t in rng.random(num_points_circle)]).T +\
    rng.multivariate_normal(mean, cov, num_points_circle).T # points as columns
disk = np.sqrt(rng.random(num_points_disk)) * np.array([[np.sin(2*np.pi*t), np.cos(2*np.pi*t)] for t in rng.random(num_points_disk)]).T
points = np.concatenate((circle, disk), axis=1)
plt.scatter(circle[0, :], circle[1, :], s=10)
plt.scatter(disk[0, :], disk[1, :], s=10)
plt.gca().set_aspect('equal')
colours = [0]*num_points_circle + [1]*num_points_disk

We compute the chromatic alpha complex \(K\) of the point cloud. \(K\) has far fewer simplices than either the Cech or Vietoris-Rips complex, which each have \(\displaystyle \binom{200}{2} = 19900\) edges and \(\displaystyle \binom{200}{3} = 1313400\) 2-simplices.

K = ch.chromatic.alpha(points, colours)
print(f'{len(K.simplices[1])} 1-simplices')
print(f'{len(K.simplices[2])} 2-simplices')

954 1-simplices
1312 2-simplices

We can visualise the filtration at any given time. For example, at time 0.2 we have:

ch.plotting.draw_filtration(K, points, 0.2)

(<Figure size 640x480 with 1 Axes>, <Axes: >)

We can compute the six-pack of persistent homology diagrams of the chromatic alpha complex, the chromatic Delaunay-Cech complex, and the chromatic Delaunay-Rips complex.

with warnings.catch_warnings():
    warnings.simplefilter('ignore')
    truncation = 1.0
    dgms_alpha = ch.sixpack.compute(points, colours, dom=[0,], method="chromatic alpha")
    fig_alpha, ax_alpha = ch.plotting.plot_sixpack(dgms_alpha, truncation)
    fig_alpha.suptitle('Chromatic Alpha')
    dgms_delcech = ch.sixpack.compute(points, colours, dom=[0,], method="chromatic delcech")
    fig_delcech, ax_delcech = ch.plotting.plot_sixpack(dgms_delcech, truncation)
    fig_delcech.suptitle('Chromatic Delaunay-Cech')
    dgms_delrips = ch.sixpack.compute(points, colours, dom=[0,], method="chromatic delrips")
    fig_delrips, ax_delrips = ch.plotting.plot_sixpack(dgms_delrips, truncation)
    fig_delrips.suptitle('Chromatic Delaunay-Rips')

Decomposed g
Decomposed cok
Decomposed rel
Decomposed im
Decomposed ker
Decomposed f
Decomposed rel
Decomposed f
Decomposed g
Decomposed im
Decomposed cok
Decomposed ker
Decomposed rel
Decomposed f
Decomposed g
Decomposed cok
Decomposed im
Decomposed ker

Each six-pack is an object having attributes ker, cok, im, dom, cod, rel, entrance_times, and dimensions. Each of the first six of those attributes is an object that has sets of paired and unpaired simplices, while the last two attributes indicate the filtration time of those simplices. For example-

dgms_alpha.ker.paired

{(228, 235),
 (252, 438),
 (262, 372),
 (280, 317),
 (282, 338),
 (292, 383),
 (331, 416),
 (333, 581),
 (345, 376),
 (348, 528),
 (352, 540),
 (358, 409),
 (359, 412),
 (361, 456),
 (387, 587),
 (413, 428),
 (421, 491),
 (429, 502),
 (436, 658),
 (454, 471),
 (487, 531),
 (506, 596),
 (544, 646),
 (584, 883),
 (606, 626),
 (1187, 1200),
 (1248, 1252),
 (1267, 1268),
 (1329, 1331),
 (1360, 1458),
 (1366, 1408),
 (1393, 1403),
 (1445, 1453),
 (1447, 1495),
 (1454, 1457),
 (1461, 1494),
 (1464, 1470),
 (1486, 1499),
 (1505, 1516),
 (1510, 1514),
 (1565, 1600),
 (1592, 1671),
 (1648, 1650),
 (1692, 1693),
 (1697, 1701),
 (1995, 1999),
 (2013, 2433)}

We can save the diagrams to file or load a diagram from file:

with h5py.File('test.h5', 'w') as f:
    ch.sixpack.save_diagrams(dgms_alpha, f)

with h5py.File('test.h5', 'r') as f:
    dgms = ch.sixpack.load_diagrams(f)

We can visualise the 2-skeleton of the filtration for points in 2D:

animation = ch.plotting.animate_filtration(
    K, points, filtration_times=np.linspace(0, 1.0, 45), animation_length=5)
animation

Once Loop Reflect